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Single-Mode Projection Filters for Modal Parameter 
Identification for Flexible Structures 




by 

Jen-Kuang Huang* and Chung-Wen Chen^ 


Abstract 

Single-mode projection filters are developed for eigensystem parameter 
identification from both analytical results and test data. Explicit 
formulations of these projection filters are derived using the orthogonal 
matrices of the controllability and observability matrices in the general 
sense. A global minimum optimization algorithm is applied to update the 
filter parameters by using the interval analysis method. The updated modal 

parameters represent the characteristics of the test data. For illustration 
of this new approach, a numerical simulation for the MAST beam structure is 
shown by using a one-dimensional global optimization algorithm to identify 
modal frequencies and damping. Another numerical simulation of a ten-mode 
structure is also presented by using a two-dimensional global optimization 
algorithm to illustrate the feasibility of the new method. The projection 
filters are practical for parallel processing implementation. 


^Assistant Professor, Dept, of Mechanical Engineering and Mechanics, Old 
Dominion University, Norfolk, Virginia 23529-0247. 

“Research Assistant, Department of Mechanical Engineering and Mechanics, Old 
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Chapter 1 
INTRODUCTION 

The problems of deriving control algorithms and state estimators for 
maneuvering flexible structures have been investigated by many researchers in 
recent years. The control design demands an accurate model of the system 
dynamics which will adequately describe the dynamic behavior of the system. 
System identification methods use experimental measurements to estimate 
dynamic properties such as natural frequencies, damping factors, mode shapes 
and modal masses which are referred to as modal parameters. Several different 
time-domain and frequency-domain methods are possible for the identification 
of structures. Various techniques may share the same mathematical foundation 
via system realization theory^. However, most techniques do not account 
explicitly for the factors which affect the actual performance in practice 
significantly. These factors include nonlinearities, local modes, and system 
and measurement noise. In order to achieve the final purpose of identifica- 
tion, i.e., control of flexible structures, an on-line estimation technique 
needs to be used. 

For linear time- i nvari ant systems, optimal model-reduction and state 
estimation has been developed via optimal projection equations based on 
modified Riccati and Lyapunov equations^. Other filtering approaches in both 
time and frequency domains are easy to implement and effective in rejecting 
uncorrelated measurement noise from simulated data^. However, previous time- 
domain filters usually involve an unacceptable computational burden for multi- 
mode systems such as large flexible structures. 

In this paper, simple projection filters are developed using system 
realization theory which has been considerably used in deriving system 
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identification methods*. The development of projection filters was initially 
motivated by the need of an on-line technique for state estimation of linear 
dynamical systems. Because of modeling errors due to system uncertainties, 
the projection filters are naturally required to be verified and updated from 
measurement data. The main objective of this paper is to present a novel 
approach to update the projection filters which in turn yields the modal 
parameter identification for linear dynamical systems. 

Each projection filter is formulated with a single mode only and its 
explicit expression can be derived using the orthogonal matrices of the 
controllability and observability matrices in the general sensed The modal 
parameters (including modal frequency, modal damping and mode shapes) required 
for formulating the projection filters are initially obtained from an 
analytical model in modal space. The experimental data are then passed 
through the projection filters to determine whether there is a discrepancy 
between the analytical model and the experimental testing. Each projection 
filter is then updated by varying the corresponding modal parameters to 
minimize a cost function defined by the norm of a specified error matrix. 
Both one and two-dimensional global minimum optimization algorithms® - ® are 
applied for the filter update by using the interval analysis method which 
guarantees finding the smallest value of a cost function throughout a 
specified closed region of modal parameters. The updated projection filters 
thus produce the modal parameter identification for the system. 

Finally, two numerical examples, one for the MAST truss beam structure 9 
and another for a ten-mode structure are given to illustrate this new method. 
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Chapter 2 

PROJECTION FILTERS FORMULATIONS 

The projection filters are developed using system realization theory. A 
finite-dimensional, linear, time-invariant dynamic system can be represented 
by the state-vari abl e equations in discrete-time form: 

q ( k+1 ) = A q(k) + R u(k) (1) 

y ( k ) = C q(k) (2) 

where q is an n-dimensional state vector, u is an m-dimensional control or 
input vector, and y is a p-dimensional measurement or output vector. The 
integer k is the sample indicator. For flexible structures, the state 
transition matrix A is a representation of mass, stiffness, and damping 
properties. The control influence matrix B characterizes the locations and 
type of input control vector u. The measurement influence matrix C describes 
the relationship between the state vector q and the output measurement vector 
y, and characterizes the mode shapes of the system. 

For the state-vari abl e Eqs. (1) and (2) with free pulse response, the 
time domain description is given by the function known as the Markov parameter 

Y(k) = CA k_1 B (3) 

or in the case of initial state response (zero input) 

Y(k) = CA k q ( 0 ) 

where q(0) represents the initial conditions of the state vector and k is an 
integer. The functions Y(k) can be obtained from the measured data and used 
to form the (r+1) by (s+1) block data matrix (generalized Hankel matrix) 
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Y(k) 

Y(k+t 1 ) ... 


H(k-l) = 

Y(j 1 +k) 

• 

• 

Y(j 1 +k+t 1 ) 

• 

• 

....Y(j 1 +k+t s ) 

• 

• 


• 

V(j r +k) 

• 

Y( j p +k+t 1 ) 

• 

Y( j r +k+t s ) 


(4) 


where j 1 -(i=l,...,r) and (i=l, . . . ,s) are arbitrary integers. For the system 
with initial state response measurements, simply replace H(k-l) by H(k). 

From Eqs. (3) and (4), it can be shown that 


H(k) = V f A K W $ ; V p = 


and 


CA 


J r 
CA r 


W $ = [B, A 1 B,..., A S B] 


(5) 


where V p and W s are generalized observability and controllability matrices 
respectively. The dimensions of V r and W $ are (r+l)p x n and n x m(s+l) 
respectively. Now observe that 


H ( 0 ) = V p W $ (6) 

If A is nonsingular, (r+1) p > n and m(s+l) > n, we can derive 

vf H ( 0 ) wf = I with Vf V n = W c wf = I (7) 

r sn rrssn v ' 

R R 

where V* and kr are the orthogonal matrices of V p and W s , respectively, 
and I n is an identity matrix of order n. Now, instead of using the matrices 

R R 

V* and W s for an n-dimensional multi-mode system, we develop simpler forms of 
V # and which represent the orthogonal matrices of the respective 

generalized observability and controllability matrices derived from a single- 
mode model only. Note that and are rectangular matrices with dimensions 
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2x(r+l)p and m(s+l)x2 respectively. The general explicit expressions of V # 

and W # will be derived later. The matrices V # and W # , which are formulated 
only for the specific modes of interest from the analytical results, will be 
used as the left and right projection filters respecti vely. The Hankel matrix 
H(0), which is formed from measured data, will then pass through the 
projection filters to identify the modal parameters which represent the 
characteristics of the measured data. If the modes of the measured data are 
uncoupled and distinct and the projection filters have the same modal 
characteristics as the measured data have, then from Eq. (7) we have 

V # H(0) W # * I 2 (8) 

where I? is a 2x2 identity matrix. The approximate equality in Eq. (8) is 

caused by the finite sampling time T used for these digital filters. If T is 
equal to zero for a finite data length, we have (see Appendix I for proof) 

V # H ( 0 ) W # = I 2 (9) 

On the otherhand, if Eq. (8) does not hold, it indicates that the modal 
parameters of the projection filters are different from those of the measured 
data. The projection filters should be tuned in order to match the modes of 
the measured data. The algorithm for the filter update is developed in the 
next section. 

Now, the explicit expressions of the projection filters V # and can be 
derived as follows. A single mode, continuous-time, linear, time-invariant 
dynamic system has the state-variabl e equations in modal space 


x = A x + B u 


( 10 ) 
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with 


y = C x 



(ID 

( 12 ) 


where oj is the modal damped frequency (the imaginary part of the eigenvalue) 

and -a is the damping (the real part of the eigenvalue). If we denote the 

w n and £ as the natural frequency and damping ratio respectively, then we 

2 1/2 

have u = u n (1 - C ) and a = C w n * The corresponding SISO discrete- 

time system can be represented by Eqs. (1) and (2) with 


and 


e _a ^ cosojT e” a ^ sin^T 

-e" aT sinajT e' aT cosaJ 



C - [ Cj 



(13) 


(14) 


where bj , b£, c^, C 2 are scalars. From Eq. (5) with j 0 = t 0 = 0 


CA 


Ji 

V = CA = [c 2 V x (r) + Cj V 2 (r), - Cj Vj(r) + c 2 V 2 (r)] 


(15) 


CA 


t t. t 

W = [A 0 B, A 1 B,..., A s B] 


-b 2 vj(s) + b x vj(s) 
bj vj(s) + b 2 vj(s) 


(16) 


where 


v[(r) = [..., V^, ...], vj(r) = [..., V 2i , ...] 


(17) 


- LctT -j.-oT 

V xi = -e sin(j.a)T), V 2i = e cos(j.coT) (18) 
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with i~0, 1, 2,*«#, r 



vj(s) = [..., v 3k , ...], vj(s) = [..., v 4k , ...] 

(19) 


-t k aT -t k <jT 

V 3k = -e sin(t k wT), V 4k = e cos(t k u>T) 

(20) 

with 

k = 0, 1, 2,... ,s 


Assume we choose j^, t k as follows 



J r - J r _i = J'i i = 0, 1, 2 ,..., integer [£] 

(21) 


t $ - t $ _ k = t k k = 0, 1, 2 ,..., integer [-|] 

(22) 

Then, the 

projection filters, V # and W*, have the following 

expl ici t 

expressions 

(see Appendix II for proof): 



( c 2 vj(r) + Cj vj (r))/(cf + c\) 

/ = 

(- c l vj(r) + c 2 vj (r))/(c \ + c\) 

* ' 

(23) 

w* - [(-b 2 

vj(s) + b x vJ(s)) T , (bj vj(s) + b 2 vJ(s)) T ]/(b^ + b\) 

(24) 


where 


v'( p ) 


[ 


• • • 
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]. 




• • • 


] 



J,aT 

e (- si n( j i wT) (l/X^r) + l/\ ? ( r) ) - sin ( (j p -J 1 )coT) 


(25) 
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(1/X 1 (r) -l/X 2 (r))/2 


= e Jl (cos(j.a)T) (l/Xj (r) + 1/X 2 ( r) ) + cos( ( j p -j i )u>T) 


(1/X 1 (r) - l/\ 2 (r))/2 


(26) 


with i = 0,1,2,... ,r 


and 


v:(s) = [..., VI, ...], Vl(s) = [..., V/, ...] 


(27) 


\l* = e k (- sin(t k ajT) (l/x 3 (s ) + l/x 4 (s ) ) - sin ((t g - t k ) w T) 
c k 


(1/X 3 (s) - 1/\ 4 (s)))/2 


Vj = e k (cos(t k u>T) (l/x 3 (s) + l/x 4 (s)) + cos((t s -t k )a)T) 
a k 


( 1 /X 3 ( s ) - 1/x 4 (s)))/2 (28) 


with k = 0, 1, 2 s 

X^(r) and x 2 (r) are the eigenvalues of V T V, and x 3 (s) and X 4 (s) are the 

eigenvalues of WW^ when = c 2 = 1 and b 2 = Cj = 0. X^(r), X 2 (r), X 3 (s) and 
X 4 (s) can be derived as follows 


X 1 (r) = 


r l+m+Y(m) 
(m+Y (m) 


if r is even 
if r is odd 


(29) 
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X 2 ( r ) = m- Y ( m ) 


(30) 


with 
m = { 


£ 

2 

r+1 


if r is even 
if r is odd 


m- 1 

Y(m) = e cos( (j -2j.) u T) 
i=0 r 1 


x (si = / 1+n+z ( n ) 

x 3 [s) W(n) 

\ 4 (s) = n-Z(n) 


if s is even 
if s is odd 


(31) 

(32) 

(33) 


with 


n = { 


s 

1 


if s is even 


s+1 ’ 

51 if s is odd 


n-1 

Z(n) = e cos( (t_-2t. )uT) 
k=0 s K 


(34) 


Note that V # and are rectangular matrices with dimensions 2x(r+l) and 
(s+l)x2 respectively. 

The corresponding MIMO single-mode subsystem can be represented by Eqs. 
(1) and (2) with the A matrix shown in Eq. (13) and 


B = 


... b 
... b. 


If *•' 


C = 


c gi c g2 


with 


2f 

f = 1, 2, ...,m and g = 1, 2, ..., p. 


(35) 


Then the observability and controllability matrices V and W can be derived as 
fol 1 ows 


V = 


CA 0 




j l 
CA 1 

= 

• 

• 

• 

V i 

W = [A 

*‘j r 

CA r J 


• 

• 

• 

II 


t t. t 

\ A 1 B, .... A ' 

[••••» ••••] 


(36) 


10 


where 


V i = 


C 12 V li + C 11 V 2i 
C 22 V li + C 21 V 2i 


C p2 V li + C pl V 2i 


- c n \i u + c 12 v 2i 

" C 21 V li + C 22 V 2i 


- C pl V H + C P 2 V 2i 


W, = 


-b 

b 


21 

V 3k + 

b ll 

V 4k’ 

" b 22 V 3k + 

b 12 

V 4k 

11 

V 3k + 

b 21 

V 4k’ 

b 12 V 3k + 

b 22 

V 4k 


5 2m V 3k + b lm V 4k 
5 lm V 3k + b 2m V 4k 


with 


i = 0 , 1 y . . . , r , k - 0, 1 9 •••9 s 


V 1 -j , v 2i , V 3k and V 4k are shown in Eqs. (18) and (20). The corresponding 
orthogonal matrices, V # and are 


where 


V # = [. . . , vf , . . .1, = 


*# 

W k 


(37) 


A- 


(c 12 V ai + C 11 V bi^ pc ll + pc 12^ * * * * ’ ^ c p2 V ai + c pl V bi^ pc pl + pc p2^ 
(~ c ll V ai + c 12 V bi^ pc ll + pc 12^’ * * * ’^" c pl V ai + c p2 V bi^ pc pl + pc p2^ 


wr = 


( " b 21 V ck + b ll V dk )/ ( mb |i + mb n) 
^ _b 22 V ck + b 12 V dk )/(mb 22 + mb 12* 


( - b 2m V Jk + b lm V dk )/ < mb 2n. + (b lm V ck + b 2m V dk )/(mb 2m + mb lV 


^ b ll V ck + b 21 V dk^ mb 21 + mb ll^ 
(b 12 V ck + b 22 V dk^ mb 22 + mb 12^ 
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V^. , V^., v£ k and V^ k are shown in Eqs. (26) and (28). Note that and W # 
are matrices with dimensions 2 x(r+l) p and (s+l)mx2 respectively for the MIMO 
case. 

In order to update the projection filters to identify the modal 
parameters from measured data within a specific range of accuracy, a cost 
function is formed as follows: 

J = j U T U (38) 

where 

U T = [E u , E 12 , E 21 , E 22 ] 

= V # H(0) W # - I 2 

From Eq. (8), the cost function J would go to a minimum value (ideally, zero) 
when the projection filters have the same modal characteristics as the 
measured data have. On the other hand, for the specified region of system 

parameters of interest, we may update those system parameters of the 

projection filters within the specified region so that the cost function is 
globally minimized. This global minimum of the cost function in the specified 
region will provide the "best" estimates for the modal parameters which 
represent the characteristics of the measured data. Although the cost 

function may be corrupted by system or measurement noise, the system 
parameters corresponding to the global minimum are expected to be quite 

insensitive to noise that is not correlated with the parameters being 
estimated. 


and 


'll 

: 12 


-12 

E 22 
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Chapter 3 

GLOBAL OPTIMIZATION USING INTERVAL ANALYSIS 

3.1 One-Dimensional Global Optimization 

The method for computing a global minimum within a specified region of 
system parameters is based on the algorithm developed by Hansen 5- ®. This 
algorithm can deal with problems in the multi-variable case with inequality 
constraints 5 . In this section, only a single variable 5 (either modal 

frequency or modal damping) is considered. Two-Dimensional case will be given 
in the next section. The global optimization algorithm basically uses a 

Q 

Newton method 0 in conjunction with the interval analysis to solve a system of 
nonlinear equations. The term "global minimum" used herein refers to the 
smallest value of the cost function J throughout a closed interval of a system 
parameter. Because of the interval analysis, the computational procedure of 
this algorithm requires explicit expressions for the first derivative (J‘) and 
the second derivative (J' 1 ) of the cost function J, Eq. (38). This can be 
derived easily by using the explicit expressions for the modal filters shown 
in Eqs. (23) and (24). The algorithm developed by Hansen 5 has been slightly 
modified. Instead of using three lists (L 0 , L^ and L 2 ), only two lists are 
applied. The algorithm is summarized as follows: 

Initial Step : The algorithm starts with an initial interval X 0 . This 

interval is equally subdivided into subintervals which are stored in a list 
Lq. A list L^ (initially empty) consists of intervals for which the width is 
smaller than a specified value w^ and the corresponding width of J is smaller 
than a specified value w^. Let x denote a feasible approximation to the 
global minimum. If the feasible point is not given, the upper limit of the 
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cost function is set to 3 = ® with x indefinite. Let [ j L , j R ]» 

[j^, j R ] and [ j * , j p ‘ ] denote the interval resulting from evaluating J, J' 
and J' 1 2 3 4 5 in interval arithmetic using the argument X, respectively. That is, 

J(X) = [j L j R ], J'(X) = [j', y R l, J"(X) = [j“, jj 1 ] (39) 

and X = [x L , x R ] 

Then, use the interval analysis to find the corresponding J, J' and J 11 for 

all the subintervals in Lq. 

Main Steps : 

1. If the list Lq is empty, go to step 11. Otherwise, find the subinterval 

X in Lq for which the left endpoint of J(X), i.e. j^, is smallest. 

2. If x e X, set x = x. Otherwise, set x = m(X) = midpoint of X. If 

j. > D, the cost of any point inside the intervals in Lq exceeds the 

upper limit D. Then Lq is set empty and go to step 11. 

3. Concavity Check 

If j R ‘ <0, J is concave in X and cannot have a minimum in the interior 
of X. Then, X is deleted and go to step 1. 

4. Monotonicity Check 

If j R <0 or > 0, the gradient of J is strictly positive or 

strictly negative over X. Then, X is deleted and go to step 1. 

5. Gaussian Elimination 

Denote E = 0 - J(x), a = [J'(x)]^ + 2Ej^' 

If j^' > 0 and a> 0, it implies that J(y)>3 for any y e X. Then, X 

is deleted and go to step 1. Note that this is true only for 1 > 0, 

which is not indicated in Ref. 5. 
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6. Interval Newton Method^ 

If j|' > 0, denote S' = x - J'(x)/J''(X) and S = intersection of S' 

L I 

and X. Otherwise, denote S = U S 2 
Here and S 2 are defined as follows: 


Denote c = 
and d = 

If J'(x)>0 


If J ' ( x ) <0 


x - 

X - 


J‘ (x)/j^- 

[x, , d] 

= I L 
1 empty 

[c, x R ] 
^ empty 


- { 


[x L , c] 
empty 


= { 


[d, x R l 
empty 


for j^' * 0. 


"Jr' 

* 

0 . 



when 

Jr 

'>0 

and 

d>x^ 

when 

Jr 

•=o 

or 

d<x^ 

when 

J L 

'<0 

and 

C<X R 

when 

Jl 

'=0 

or 

ox R 

when 

K 

'<0 

and 

C>X L 

when 

K 

'=0 

or 

c<x^ 

when 

j ‘r 

•>o 

and 

d <X R 

when 

Jr 

'=0 

or 

d>X R 


7. If S is empty, go to step 1. If S=X, then split X in half. 

A 

8. For each new generated subinterval, X = or S, repeat steps 3 

and 4. 

9. Update 0 

A 

For each new subinterval X, denote 

w[X] = width of X, x = mid [X] = midpoint of X 

X = [x L , x R ] and J(X) = [j,_, j R ] 

A A 

If J(x)<D, simply replace 0 by J(x) or conduct a line search to 


reduce J as follows: 
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a. If J'(x)>0, denote x^ = x L « Otherwise, denote x^ = x R . Set 

A A 

x o = x * 

AAA A A A 

b. Denote x ? = (x Q + x^/2. If J(x 2 ) > max [J(x Q ), J(x 1 )], go to 
step e. 

A A A A A A 

c. If J(Xq) <J(x 1 ), replace x^ by x^. Otherwise replace x Q by x^. 

• A A 

d. If | - x Q | > w[X], go to step b, 

A A A 

e. Set 3 = min[J(x), J(x Q ), J(x^)] and set x to the corresponding 

arugument of 0. 

10. Store new intervals 

A A A 

For each new interval X, if > 0, delete X. 

AAA A A 

If w[X] < w 1 and j R - j' L < w^, store X in Lj. Otherwise, store X 
in Lg. Go to step 1. 

11. If the list Lj is empty, go to step 13. Otherwise, delete subintervals X 
for which j|_ > 3 where J(X) = [j|_, j R l. 

12. If the list I.J is empty, go to step 13. Otherwise, the midpoints of each 
interval remaining in are used as the global minima. Note that there 
may exist multiple global minima. 

13. Because Lj is empty, the global minimum is located on one of the two 
boundaries of the initial interval X Q , which corresponds to a smaller J. 

3.2 Two-Dimensional Global Optimization 

3.2.1 The Statement of the Problem 

Given a scalar objective or cost function f(x) of two independent 
variables, find a point x* within a specified rectangular region X^, which 
is called "box" in the sequel, X^ R^, such that 
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f(x*) < f(x) (40) 

for any x e X ^ . 

There is no constraint except that x must be within a rectangular region 
with sides parallel to the coordinate axes. Here, we assume the objective 
function is twice continuously differentiable. 

The algorithm presented here is composed of four main separate parts. 
Let X be a sub-box of X^) and is the box under treatment. One part uses an 
interval version of Newton's method to find a sub-box or some sub-boxes of X 
such that any stationary point within X will still be contained in these sub- 
boxes. By this method, we can shrink the box under treatment to a smaller box 
or several smaller boxes without losing any possible extreme point. A second 
part eliminates points of X where f is greater than the smallest currently 
known value. This is because only the global minimum is of interest. Any 
point which is obviously not the minimum could be deleted. A third part tests 
the monotonicity of the box. If f in the box X increases or decreases 
montonical ly, this box can be eliminated, except the box contains the boundary 
points of the initial box x(0). a fourth part checks the convexity of the 
sub-box X. If f is not convex anywhere in X, there cannot be a stationary 
minimum of f in X. Therefore, the whole box can be deleted. Note that if the 
global minimum point x* occurs on the boundary of X^), it is not necessary to 
be a stationary point. So when we use the fourth part of the algorithm, we 
must retain all the boundary points. 

The following several sections will describe these four parts in more 


detail 
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3.2.2 The Interval Newton's Method 

This method seeks the zeros of the gradient function g of the objective 
function f and hence the stationary points of f. In ordinary functions, by 
using Newton's method, we can locate the stationary points of f. In interval 
analysis, by using the interval version of Newton's method, we can find a new 
box or several boxes N(X), which are smaller than the previous box X, such 
that any point in X not in N(X) cannot be a zero of g and can be discarded, 
unless it is a boundary point of X^). The method is described as follows. 

f(x) is a scalar objective function, x is a column vector, and g(x) is 
the gradient function of f. Note that g(x) is also a column vector. We want 
to solve the equation 

g ( x ) = o (4i) 

Let y be an approximate solution to Eq. (41) and J(x) denote the Jacobian 
matrix. Then using Taylor's theorem to expand g(x) about y, we can get the 
following equation 


g( x) = g(y) + J(x) (x-y) = 0. (42) 

where x is a point between x and y. If x, y e X then x e X. We will obtain 
a box or several boxes as its solution if substitute real point 
x by a box X and solve Eq. (42). 

To solve Eq. (42) with x replaced by X, first we find a real matrix J c 
which is the center of the interval Jacobian matrix J(X). That is, each 

element of J c is the center of the correspondent interval element of the 
interval Jacobian matrix. Let B be an inverse of J c , then on 
premultiplication of Eq. (42) by B, we get 


Bg(y) + BJ(X) (x - y) = 0 


(43) 


Any point in X solving Eq. (41) should be contained in the interval solution 
of Eq. (43), which could be solved by the following method. Write 

B J ( X ) = L + D + U (44) 

where L, D and U are the lower triangular, diagonal and upper triangular part 
of BJ(X) respectively. The interval matrix 

D" 1 = diag [1/D n , 1/D 22 ] (45) 

is the inverse interval matrix of D. Then the solution of Eq. (43) is 

X' = y - D' 1 [Bg (y ) + L(X* - y) + U(X - y)] (46) 

Though X 1 also appears in the right hand side of Eq. (46), its elements can be 
solved recursively, because it is multiplied with a lower triangular matrix 
L. If 0 i and 0 l X' is a single box, otherwise X' could be several 

boxes. After intersecting with the previous box X, the answer could be the 
empty set, one box which is the same as X or a reduced box or several reduced 
boxes. 


3.2.3 Bounding ? 

As we proceed with the algorithm, we should evaluate f(x) at a number of 
points of X(°). Let T denote the currently smallest value of f found so far, 
and let X be any sub-box of x(^). We can delete any point x of X where f(x) 
> T. The method we used is described as follows. 

Expanding f about a point x e X, one obtains 


f(y) = f ( x) + (y - x) T g(z) 


(47) 
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where z is a point between x and y, g is the gradient function of f. If 
x e X, y e X, then z e X. Now we substitute z by the box X. If 


f(x) + (y - x) T g ( X ) > T 


(48) 


then f(y) > ?. We want to retain points where Eq. (48) is not satisfied. 
This is the same as to say we want to retain those points where the following 
inequality is satisfied. 


f(x) + (y - x) g(X) < T 


(49) 


Let 


then 


f - f ( x ) = e 


(y - x) g(X) - e < 0 


i .e. 


Let 


then 


C(y 1 - * 1 )> (y 2 - x 2 )] 


y\ - x i = h 


g x (x) 

g 9 (x) 


- e < 0 


y 2 - x 2 = y 2 


y! gi (X) + y 2 g 2 (x) - e < 


0 


(50) 


If we expand the function about the center of X, and first use X 2 of X(X^, X 2 ) 
as y 2 , then Eq. (50) has the general form of a linear equation with interval 
coefficient 

A + Bt < 0 (51) 
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Here, = t, A = X 2 g 2 (X) - e, B = gj (X). Let A = [a lt a 2 ], B = [bj, b 2 ], T 
= {t: at + b < 0, aeA, beB}, then the general solution T to Eq. (51) is as 
fol 1 ows 


[-a 1 /b 2 , -] 
[■ /b^ * 00 3 


T = < 


[- 00 » ® ] 
[-«, -a 1 /b 2 ] 

c~ 00 9 ~ ^ /bj ] 

[- 00 * - a^ /b 2 3 


U [-aj/bj, «] 


empty set 


if a^ < 0, 
if a^ > 0, 
if a^ < 0, 
if a^ > 0, 
if a 1 < 0, 
if a 1 > 0, 
if a^ > 0, 


b 2 < 0 

b^ <0, b 2 < 0 
b^ < 0 < b 2 
b^ < 0 < b 2 
b l > 0 

bj > 0, b ? > 0 


(52) 


Note the solution of Eq. (50) may be empty set, one interval or two 
intervals. Though the interval may be semi-infinite, after intersecting with 
the previous box X, the result is finite. After obtaining Y^, we treat Y 2 as 
a variable and solve it similarly. 


3.2.4 Check for Monotonicity 

Let X be any interior sub-box of X^), which means X contains no boundary 
points of X<°>. Note if any optimal point should occur in the interior it 
must be a stationary point. Therefore, if we find f increases or decreases 
monotonical ly throughout X with respect to each variable x^, we can delete the 
whole box X, since X does not contain any stationary point. To make use of 
this characteristic, we first evaluate the interval gradient function g of f 
and check each component. Denote [d j , Wj] for the interval value of the j-th 
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component g j ( X ) of g, (j = 1, 2). It is obvious that dj < g j ( x ) < Wj for any 
x e X. Denote the i-th component of X by X^ = [x^, x^], if dj = 0 and wj > 
0, in the x,- direction f(x) is smallest for x. = x^ (in X). So all of X 
except the boundary line of which Xj is its component can be deleted. If dj > 
0 we can delete all of X, except if Xj is a boundary point of X^). in this 
case, we save this boundary point. Similar results can be obtained if w^ = 0 
and d i < 0 or if w i < 0. We check g^X) (i = 1, 2) one by one, delete points 
as many as possible and save those that cannot be deleted. By this method, we 
may delete the whole box, reduce the box to a degenerated box, a line segment, 
or reduce the box to a point. In the latter case, we evaluate f at this 
point. If f > ?, this point can be deleted. If f < ?, we replace ? by f and 
store it. 


3.2.5 Test for Convexity 

If f is not convex in a given box, none of its stationary points can be a 
minimum. Because we are only interested in the global minimum, if f in a sub- 
box X of X<°> is not convex, we just delete it before doing any other 
treatment. But any boundary point of the initial box X^) should be 

retained, for if the minimum occurs in the boundary it is not necessarily a 
stationary point. 


To check the convexity of a function, we need to check its second 
derivative. For a scalar function f of two variables, the second derivative 
is a 2 x 2 Hessian matrix 


H 


a 2 f 

aXj 

a 2 f 


ax 2 aXj 


a 2 f 


aXj ax 2 



(53) 
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For a real function, the stationary point is a minimum if the Hessian matrix 
evaluated at this point is positive definite. For interval function, it is 
not easy to check if the Hessian matrix is positive definite. But we can use 
a quick way to rule out those points which are definitely not qualified. This 
can be done by checking the diagonal elements of the interval Hessian 
matrix. For a positive definite square matrix, each diagonal element must be 
greater than zero. Taking advantage of this fact, we first evaluate the 
interval value of the diagonal elements. Denote the i-th diagonal element by 
d^-fx) and its corresponding interval value of the i-th diagonal element 

by CnV ^ , D*?.](i = 1» 2), if D^ < 0 then d^(x) < 0 for any x e X, the Hessian 

matrix can not be positive definite and we can eliminate X, except those 
boundary points of X^). 

3.2.6 Boundary Points 

If the global minimum is a stationary point, then we are free to delete 
all the boundary points. Otherwise, we should be careful not to delete 
boundary points which might be the global minima within the specified 
region. The interval Newton method and the procedure which checks the 
convexity cannot delete any boundary point of X^ 0 ^. In these two procedures, 
if we find that the treated sub-box contains boundary points, we just bypass 
the procedures. However, the procedures of checking montonicity and 

bounding ? as described above can eliminate some boundary points. When 
checking monotonicity, if we find g^ (Xj, X 2 ) > 0 (i = 1, 2) then we can 

delete all points of X except if the points at the left boundary line of X 
with respect to are also the boundary points of X^) # In this case, we 

save this boundary line and delete all the other points including some 

boundary points of X^ 0 ) which might exist. If g i (Xj, X 2 ) < 0, similar 
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treatment can be applied. When bounding ?, we can delete any point (including 
boundary point) x if f(x) > ?. 

We can also separate the boundary points from the interior points of 
in the beginning of the process, and treat the interior as a box wherein the 

global minimum must be a stationary point. Those separated boundary points 

are processed later. This procedure is more straightforward but will be 
impractical when the dimension of is getting larger because it will 

generate too many boxes. 

3.2.7 The Lists of Boxes 

Our algorithm begins with an initial box X^). By processing X^^ 

through the above procedures for one run, we will either get an unchanged box, 
one smaller box or several sub-boxes. For the unchanged box, we equally 
divide it by two and save them in a list denoted by L^. For the smaller box, 
we check how much it was reduced. A criterion was set to make a judgement. 
If it is only slightly reduced, we divide it by two along its longest 
dimension and save them in L^. If it has been reduced a lot, then just save 
it in Lj. For the case of several sub-boxes, we save them all in for 
further treatment. At the beginning, the number of boxes saved in will 
increase rapidly. But when the size of the box becomes smaller, we will have 
more chances to delete the whole box through the procedures described above. 
Then the number of boxes in Lj will begin to decrease. 

Before saving any box in at the end of each run, we will check whether 
the size of the box and the range of the objective function corresponding to 
the box are smaller than the specified accuracies of the solution (e^ and 
described in the next section). If it does, we save it in another list 
denoted by L ?. Eventually, the number of boxes in will become zero and 

only one or a few boxes will be left in l_ 2 . 
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3.2.8 Termination Criterions 

Two numbers, e^ and e 2 , are used to specify the accuracies of the answer 
we want. By which we mean the width of each dimension of the solution box 
must be smaller than e^ and the interval of its objective function must be 
smaller than e 2 . If we denote the width of the i-th element by w(X^) and 
the interval of the objective function by [F*-, F^], then the requirement can 
be expressed by 


max [w(X.)l < e-i (54) 

1 * 1,2 1 

F R - F L < e 2 (55) 

After l_i becomes empty, if there are a few boxes, say s, remaining in L 2 , 

we have to continue to work on them to yield the best answer. Oenote the 

boxes by X^ (i = 1, 2,...,s). For each 1 = 1, 2,...,s, evaluate f(X^)) and 

denote the result by [F*r, F R ]. Then for every x e X^, F 1 : < f(x) < F R . Now 

denote 





f = min 

F i 

s 



(56) 




1 < i < 




then 



f < f(x) 




(57) 

for x in 

any of 

the boxes X^ ) 

(i = 1, 2, . 

..,s). Since 

any 

global minimum x* 

must be 

one 

point within 

a certain 

box among 

x(i) 

, we can 

i nfer 

that f < 

f(x*). 

From the 

preceding 

section, ? 

is 

an upper 

bound 


of f(x ), i.e. ? > f(x ) for anytime, hence 

f < f* < T 


(58) 
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If e = T - f is smaller than the desired value e?, we terminate the 
algorithm and list out all the remaining boxes in I 2 as our answers. If e is 
bigger than then we find the value j for which f = F*r, remove this box 

from I 2 to Lj and re-run the main algorithm. This procedure can always 

increase f and hence decrease e. We repeat this process until e is smaller 
than ej. Some boxes in l _2 could be deleted by this procedure. Those 

remaining in L 2 are the answers. Multiple answers are possible. 

3.2.9 Steps of the Algorithm 

We now describe the steps involved in the algorithm. First, the box X^) 
is divided into several sub-boxes to reduce its size because if a box is too 
big, it will usually remain unchanged after passing through all the steps. 
Then we evaluate the value of the objective function at the center of each box 
and choose the smallest one as the initial T. The value of ? will keep on 
being reduced during the process. The subsequent steps are to be done in the 
following order except as indicated by branching: 

1. Choose the first one box remaining in l_j , call it X. If list is empty, 
go to step 8. 

2. Check for montonicity. Evaluate the gradient interval function g(X). If 

L R 

g ^ ( X ) > 0 (or < 0) and the boundary of X at x n - = x^ (or = x^ ) does not 
contain any boundary point of X^, delete X and go back to step 1. If 
g-j (X) > 0 (or < 0), replace X i = [xV, xj] by [xV, x * 1 2 3 4 :] (or [x 1 ?, x^]). 

3. Test for non-convexity. If X contains boundary points of X^), go to step 

4. Otherwise evaluate h^- (X^, X 2 ), (i = 1, 2). If either h^ or is 
strictly negative, delete X and go back to step 1. 
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4. Bounding ?. Evaluate g(X) according to Eqs. (50) and (51) and solve for 

interval Y. Intersect Y with X. If the result is empty, go to step 1. 

If the result is a single box, rename it as X. If it ends up with two 
boxes, save one in l_j, rename the other one as X and go on. 

5. Interval Newton method. Evaluate the interval Jacobian matrix and solve 
for X 1 according to Eqs. (43) to (46). If the denominator interval 
contains zero, the element Xi in X 1 becomes two boxes. Then save the gap 
between the two boxes and denote it by g^. Find the smallest box which 
contains these two boxes and rename it as X'.. Finally, intersect X 1 with 
X and rename it as X. If the result is empty, go back to step 1. 

6. Check the gap saved. Choose the largest one, say gj, and divide the box X 

along gj. If there is no gap saved, check the size of the current box 

X. If it does not satisfy the accuracy requirements e^ and & 2 > divide it 
into two boxes along its greatest dimension and go on. 

7. Evaluate f at the center of each of the boxes before storing them. 

Update T by choosing the smallest one from ? and f(x) at the center of 

each box as the new ?. If the box satisfies the accuracy requirements, 
store it in l^. Then store all the remaining boxes in and go back to 
step 1. 

8. Evaluate f(X^)) for each box X^) remaining in l^. Denote the result 
by Cf!t, F*] and find 

f = min F 1 : (59) 

‘ l<i<s 1 

If e = T - f > e^, put the box Xj which has f as its left end-point back 
to L| and go back to step 1. If e < e^, stop the process and print out 
all the boxes remaining in Lo as the answers. 
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The algorithm described above can be expressed roughly by a chart shown 
in Fig. 3.1. 


3.3 Numerical Examples 

In this section, we use two examples to illustrate the application of the 
two-dimensional global optimization algorithm. In case one, we consider the 
so-called three hump camel function 

f(x) =2xj- 1.05 xj + ^ x® - x x x 2 + x^ (60) 

which has three relative minima and two saddle points. Among them, there is 
only one global minimum. The gradient g(x) has two components 


g^(x) = 4 x^ - 4.2 + " x 2 (61) 

g 2 (x) = 2 x 2 - x 1 (62) 

The Hessian matrix of f has four components 

h u (x) = 4 - 12.6 x^ + 5 xj (63) 

h 12 (x) = -1 = h 21 (x) (64) 

h 22 (x) = 2 (65) 


The interval extension of the gradient function and the Hessian matrix 
can be derived by substituting the interval variables and X 2 for the real 
variable x^ and X 2 . The subroutines for evaluating real and interval value of 
the functions are written. Running the two-dimensional global optimization 
program by using the initial intervals X^ = [-100, 100], X 2 = [-100, 100], the 
global minimum is shown in Table 3.1 

It takes 2.8 cpu sec to obtain the answer in CDC Cyber 830 computer. 
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In case two, we consider a function called Himmelblau function 

F(x 1 , x 2 ) = - (x 2 + x 2 - ll) 2 - (x x + x 2 - 7) 2 (66) 

The isometric representation and the corresponding plot of its contours or 
level curves are shown in Fig. 3.3 and Fig. 3.4 respectively. In Fig. 3.3, 
the flat base represents a function value of -150. This function has four 
peaks and several saddle points. The global optimization problem is to find 
out its global extremes. Though the program is written for searching global 
minimum and the function here has maxima to be located, the problem can be 
easily solved by setting f(x^, X 2 ) = -F(xj, X 2 ) and find global minimum of 


f(x|, X 2 ). Now 

f(x r x 2 ) = (x 2 + x 2 - ll) 2 + (x x + x 2 - 7) 2 (67) 

The gradient g(x) has two components 

g x (x) = 4 (x 2 + x 2 - 11) x 1 + 2 ( x^ + x 2 - 7) (68) 

g 2 (x) = 2 (x 2 + x 2 - 11) + 4 (xj + x 2 - 7) x 2 (69) 

The Hessian matrix has four components 

h n (x) = 12 x 2 + 4 x 2 - 42 (70) 

h 12 (x) = 4 x 1 + 4 x 2 = h 21 (x) (71) 

h 22 (x) = 4 + 12 x 2 - 26 (72) 


We write the subroutines for calculating the value of real function and 
the corresponding interval functions and then run the program with an initial 
box of X.j = [-100, 100], X 2 = [-100, 100], The results are listed in Table 

3.2. 
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This result is exactly the same as the result found in Ref. [12]. It 
takes 2.82 CPU sec to get the answers in CDC Cyber 830 computer. 
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Chapter 4 

MODAL PARAMETER IDENTIFICATION USING PROJECTION FILTERS 

In this chapter, the global optimization algorithm using interval 
analysis will be applied to the single-mode projection filters for modal 
parameter identification. A numerical simulation of a ten-mode structure is 
first presented to demonstrate the feasibility of the projection filters as 
well as the two-dimensional global optimization algorithm. Then another 
numerical simulation of the MAST beam structure is shown by using a one- 
dimensional global optimization algorithm. The results show that the 
projection filters are feasible for modal parameter identification. 

4.1 Numerical Simulation for a 10-Mode Structure 

As described in Chap. 2, a finite-dimensional, linear, time-invariant, 
continuous time dynamic system can be expressed by the state-vari abl e 
equati ons 


x(t) 

= Ax(t) + Bu(t) 

(73) 

y(t) 

= Cx(t) 

(74) 


After appropriate linear transformation, matrix A could be expressed as 
f ol 1 ows 

A = 


■a l Wl 
-u)^ a 2 


” CT^ 

“(1)2 “0^2 


( 75 ) 


-a u) 
m m 
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which is a block diagonal matrix. Here, co^ and are respectively the damped 
modal frequency and the modal damping of the i-th mode of the system. For an 
m-mode system, the dimension of A is 2mx2m. For the discrete-time system, 
matrix A becomes 


1 


0 


* A 


A i " 


-a^T -OiT 

e cos(w.jT) , e sin^T) 

-aj -cJ.-T 

- e sin(u. T) , e cos(uj.T) 


(76) 


i = 1 , 2 , ... m 


where T is the sampling time interval. Note that each 2x2 block in the main 

diagonal corresponds to one of the system modes. Next we want to construct 

the impulse response function y according to Eq. (3). Note that 

✓ 


'1 


0 


m 


and note that 


a ; - 


-ka^T -ka-T 

e cosfkaj^T) , e sin(ko^T) 

-ka i T -k ai T 

-e sin(ku.T) , e cos(k w .T) 


(77) 


(78) 


If the system is a single-input single-output system and we denote matrix B 


and C as 
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B = [b^, ^ 2 9 •••» b m ] 


b i ' [b il, b i2 ] 


( 79 ) 


C c^, • ••» 


c i ^ c il’ c i2- 1 


(80) 


i = 1, 2, . .., m 


then 

y(k) 


C A k_1 B 


[ci , C£, • • • » c^l 


m k-1 

i?l 'i Ai b i 


* 



A k-1 

0 


b i 

A k_ 1 
A 2 

• 


b 2 

0 

•h 

• 

• k-1 

A m 

m 


• 

• 

• 

b 

m 


(81) 


m -(k-1) a H T 

E [e (c n cos( (k-l)aj^ T) - c^ sin ((k-1) u 4 T)) b.. 


'i2 


J i " il 


-(k-1) a.T 

+ e (c n sin( (k-1) w.T) - c.. 0 cos((k-l) w.T)) b. 0 ] 


i ' i2 


i ' i2- 


By assigning modal frequency and modal damping to each mode, the free 
decay impulse response function y can be constructed from Eq. (81). A series 
white noise signals are generated from a computer program and added to y to 
form a sequence of simulated output data. Using these output data, we can 
form the generalized Hankel matrix from Eq. (4). The explicit expressions of 
the single-mode filters are described by Eqs. (23) and (24). 
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Assume we have a ten-mode system with the following separate frequencies: 

3, 7, 11, 17, 23, 31, 37, 43, 53, 67. (rad/sec) 

The damping ratios are assumed to be the same for each mode. For a large 
flexible structure system, the damping coefficient is always very small. To 
simplify the problem, matrices B and C are assumed to be as follows 


B = [1, 0, 1, 0, 

i ni T 

» l » U| ix20 

(82) 

C = [0, 1, 0, 1, 

...., 0, ll lx 20 

(83) 


^ s , j .j tj i*i 1,2, ...., r 
Then the projection filters, Eqs. (23) and (24) become 


r = 


W ff = 


V a(r) 


Vj(r) 


[VJ(S) T , VJ(S) T ] 


(84) 


(85) 


(r) = Ms), \ 9 (r) = \.(s) 


V*(s) = V*(r), V*(s) = V*( r) 


Now the cost function can be expressed as follows. 


J(a, w) = + 0.5 * (E ? - ? + E^) 


( 86 ) 


where 


E n - V' H(0) Vj - 1 


Ej 2 = V* H(0) v' 
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*21 = < H <°> < 

E 22 = V' H(0) V* - 1 

Note that E n , E 12 , E 21 , E 22 are all functions of w and a. 

The first and second partial derivative of J can be derived as follows 


dJ = 9 r . 5E 

do) 8w 11 doo C 12 


12 


aE 


21 


aw 21 


(87) 


a 2 J _ 2 ^li r 

7T ‘ 2 —r E ll 

5<i) 5oo 


&E 11 2 ^ E 12 
+ — r E i2 


5o) 


a E i2 ‘ 
+ f — —) 

1 aw > 


+ (- 


a 2 E. 


aw 


ill E 

2 > l 21 


a e 2 i 

+ f— 

1 aw > 


( 88 ) 


— = 2 ^ F + 

a a a a 'll a a 


dE 12 r 8E 21 f 

fc 12 aa 


21 


(89) 


a 2 J 


a 2 E 


aa 


aa 


11 E 

T~ ni 


aE, 2 a 2 E. 

t2 h^> * 


a a' 


12 E 
T E 12 


a E-i 2 

+ f — — — ) 
aa > 


2 


a e 21 ^ e ?i 

* (— F) *2, * (-#) 


5o) 


(90) 


a 2 J = 2 

awda aw a a 


a E i i u *-i 1 

E - + 2 (-as 1 ) (-sr) * 


’ll 


5 Jj 

aa 


a 2 E 


12 


aE 


12 


aE 


12 > 


awaa 12 


^ aw ; ^ ' 


da 


a 2 E 


21 


aE, 


21 


aw3a 21 


aE, 


21 , 


+ r — LL) r — £1) 
aw ] K 1 


a a 


(91) 


with 


3E 


11 


aw 


r h(o) 

a 


a(v b > « 

+ Vf H(0) — ^ 
b aw 




5a) 


( 92 ) 
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e>E 


#T 

1? * »(V* ) 

= 2 vj H ( 0 ) — 

dco a dco 


(93) 


5E ?1 M &(< ) 

— — = 2 vf H(0) ^ — 

5 co b 


8 2 E. 


dco 


a 2 (vf) 


(94) 


= V* H(0) ' + V* H(0) 

3U) &0) 


5 2 (vj ) 

dco 2 


+ 2 


»( v f> 

dco 


H ( 0 ) 


8(<2 

5u) 


9 2 E. 2 8(V*) a(V* T ) , 8 2 (v' T ) 

r- * 2 ^T- H <°) -T5T- + 2 V a H <°> — f- 

ou) do) 

a 2 E 2 i (vj) d(vJ) T a 2 (vJ T ) 

^ = 2 -A- H ( 0 ) , D + 2 vj H ( 0 ) S— 


5o) 


50) 


(95) 

(96) 

(97) 


dco ~ w ~ dco 

The first and second partial derivatives with respect to a are all the same 
except substituting a for co. The cross partial derivative terms are 


d 2 E. 


: n av* a(v') T a 2 (v') T 

SST- = W- H <°> -T=- * U a H '°>-aJ^ 


av* a 2 (v*) T a 2 (v?) T 

a^ H(0 > ssr * < H <°> -dr 

T 

dV* 

-12 ° v a 


d 2 E i2 QV* d(V*) d 2 (vV 

a^ ■ 2 a^ H(0) + 2 u a H <°> -5ST 

# J 

- *1 _ 2 ~ y b 

5o)3a 


avf a(vj) 


:vj) , a 2 (vj) 

a» t 2 v b "«» -555T 


In numerical simulation, the following data are used 


(98) 


(99) 


(100) 


r = s = 49 


m = 10 
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T = 0.04 sec 

which means the dimension of the Hankel matrix is (49+1) x (49+1) and there 
are ten modes in the system. 

A FORTRAN program was written based on the above information. We specify 

e x = 0.001 
&2 ~ 0.01 

where e^ and e 2 are the desired accuracies of the identified system parameters 
(either frequency or damping) and of the cost function respectively. The 
noise level is the ratio of the noise standard deviation with respect to the 
maximum value of y(k), i.e. the peak of the free impulse response. Specific 
modal frequencies with different damping coefficient and noise level are shown 
in Table 4.1. Corresponding to each modal frequency, ten frequency intervals 
are given for the projection filters: 

[0.1, 5], [ 5, 10], [ID, 15], [15, 20], [20, 25] 

[25, 35], [34, 40], [40, 50], [50, 60], [60, 70] 

in rad/sec 

For each frequency interval, the range of the damping ratio is specified 
to be zero to five percent. The frequency interval along with the damping 
interval form a rectangular box. Within this box, the global minimum of the 
cost function generated by the projection filters is to be located. 

The time history of the response y(k) with 30% noise is shown in Fig. 

4.1. The white noise with 30% noise level is shown in Fig. 4.2. Figure 4.3 
is a cross section view of the cost function vs. updated frequencies as 
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damping equals to zero. Two values of measurement noise (0% and 30%) are 
shown to illustrate the effect of the noise. The result shows that the cost 
function is distorted by the noise, but the minima are not affected 
significantly. The percentage errors for the estimated modal frequency and 
damping are listed in Table 4.1 and Table 4.2 respectively. From Table 4.1, 
as the damping factor varies from 0.3% to 2%, the errors of estimated modal 
frequencies fall within 1% for the noise free case. As the noise level 
increases, the error increases proportionally and stay within 3% for 30% 
noise. For a fixed noise level, the errors increase for most of the modes as 
the damping factor increases. This is caused by the fact that the signal will 
damp out faster for higher damping factors. From Table 4.2, similar results 
are found for the modal damping errors except that the percentage errors are 
higher. For low frequency and low damping modes, the relative errors are 
higher because the contribution of the damping is comparably smaller. 

4.2 Numerical Simulation for the MAST Ream Structure 

For the multi-input multi-output system, a numerical example is given 
here based on the finite element model of the MAST truss beam structure as 
shown in Figure 4.4 (see Refs. 9 and 11 for detailed description). There are 
four actuators and 68 displacement sensors (four on each bay) distributed 
along the MAST truss beam structures’ll. There are five modal frequencies 
(see Table 4.3) and mode shapes derived from the finite element model. The 
mode shapes for the first and the second mode are shown in Figures 4.5 and 
4.6, respectively. The odd-numbered sensors measure the deflections in one 
direction of bending, whereas the even-numbered sensors measure the 
deflections in another direction of bending. Note that for these two repeated 
modes, the mode shapes are orthogonal to each other. The C matrix used for 
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the simulated data in Eq. (3) and for the projection filters in Eq. (35) is 
obtained from the five mode shapes derived from the analytical model. The B 
matrix used for this numerical simulation is arbitrarily determined from 

Y(1 ) = CB 
B = C + Y(l) 

where C + is the pseudo-inverse matrix of the C matrix and Y(l) is the first 
test data. 

The numerical example is illustrated by using the following parameters 

(see Eqs. ( 13 )- (35) ) : 

m = 4, p = 68, r = 3, s = 99, M = 5, T = 1/75 sec. 

Specific modal frequencies with different damping ratios and noise levels are 
shown in Table 4.3. The noise level is the ratio of the noise standard 

deviation with respect to the maximum value of Y(k), i.e. the peak of the free 
impulse response. For each case, the simulation starts by forming a Hankel 
matrix for this five-mode structure with a damping factor for all modes and a 
specific noise level. The simulated free impulse response data with 10% noise 
and zero damping are shown in Fig. 4.7. For each modal frequency, the 

frequency interval given for the projection filters is 0.1 to 10 Hz. With a 

fixed zero damping, the projection filters first update their frequencies by 
using the interval analysis method to find the global minimum of the cost 
function within the frequency interval. The midpoint of the final frequency 
interval (width is smaller than 0.001) is used for the first estimate of the 
modal frequency. With this estimated frequency, the projection filters then 
update their damping ratio with an initial interval from 0 to 5% by using the 
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same interval analysis method to determine the damping from the Hankel matrix. 
The midpoint of the final damping interval (width is also smaller than 0.001) 
is used for the first estimate of the damping. With this new damping, the 
whole procedure is repeated. Since the second estimates of the modal 

frequency and damping are quite similar to the first estimates, further 
estimates are prohibited. The percentage errors for the second estimates of 
the damped modal frequency and damping are then calculated for each mode and 
listed in Tables 4.3 and 4.4. The cost function J for the first mode is 

plotted in Fig. 4.8 as the projection filters updated their frequencies with a 

fixed zero damping for noise free case. The cost function is almost unchanged 

when 10% noise is added. From Table 4.3, the errors of estimated modal 

frequencies fall within 1%. For repeated modal frequencies, the projection 
filters are shown to be effective to identify those frequencies. For a fixed 

noise level, the errors increase for most of the modes as the damping factor 

increases. This is caused by the fact that the signals decay faster for 

higher damping factors. However, there are some modes for which the estimates 
of the frequency improve slightly when the measurement noise level is 
increased. For a fixed data length, if the sample time T is reduced, the 

errors in frequency decrease for the noise free case and increase for higher 
noise level. From Table 4.4, the modal damping errors are much higher as 
compared to the modal frequency errors. Because the damping ratios of the 
large flexible structure are low, they are harder to be identified. As a 

result, this numerical simulation shows that the projection filters are 
feasible for estimating the modal frequencies and damping. The global 
optimization algorithm using interval analysis is also shown to be effective. 
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Chapter 5 

CONCLUDING REMARKS 

Projection filters are derived for possible application to the state 
estimation of linear dynamical systems. An approach to update the projection 
filter through the use of measurement data is developed to identify frequency 
and damping of the system. Numerical results for the MAST truss beam 
structure demonstrate that repeated modal frequencies can be identified within 
1% error for 10% measurement noise. 

There are two characteristics of the projection filters. First, each 
projection filter is developed based on a single-mode subsystem which identi- 
fies only one modal frequency and one modal damping within a specified region. 
Second, for an n modes structure (based on the analytical model), n single- 
mode projection filters can be implemented for parallel processing to reduce 
the computational burden. Application of the projection filters to the state 
estimation for linear dynamical systems is currently under investigation. 
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APPENDIX I 


For a finite data length and a decoupled distinct M modes system as shown 
in Eqs. (37) and (38) without a measurement noise, it can be verified that 


lim V? H ( 0 ) W? = I 9 
TVO J J c 

# # th 

where V. and W. are the projection filters for j LM mode. 

J J 


(Al) 


On the other hand. 


V? V. w. W? 
J J J J 


I, 


(A2) 


where Vj and Wj are the respective generalized observability and 
controllability matrices for j th mode subsystem. 

From Eqs. (5), (6) and (A2), it is shown that 

M 

V? H(0) W? = V? V W U* = 2 V? V. W. W* 

J v J J r s j k “j j k k j 


M a u 

r 2 + k E =1 V j V k W k W j 

k^j 


(A3) 


Now, with the aids of Eqs. (15) and (23) and using the subscripts j and k to 
indicate the corresponding values for the j th and k th mode respectively, one 


obtai ns 


V # V, 
J k 


a ll a 12 
a 21 a 22 


(A4) 


a ll (c 2j c 2k V aj V lk + c 2j c lk V aj V 2k + c lj c 2k V bj V lk 


c lj c lk V bj V 2k )/(c ij + c 2j ) 


where 
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with the aids of Eqs (17), (25) and ( 29 ) - ( 3 1 ) , one derives 
u r ji (a-j-cfjJT 

V". V.. * E e J (sin(j. oj.T) s1n(j. u.T) 

a J 1K i=0 J 1 K 

r j i (a i -o k )T 

+ sin(j. w k T) sin((j r -j i ) Uj .T))/2x lj . + .E e J (s1n( J 1 oojT)si n( j. 

- sin(j. w k T) s1n((j p -j 1 )ttfjT))/2\ 2 j = 


If r is odd 

( r - 1 )/2 

X, • « ( r+1 )/2 + E cos((j -2j. ) W .T) = (r+l)/2 + Y(m) 
iJ i=0 r J 

Assume we choose t k as follows 


j i l"' ^ i » i 0, 1 , 2, . . . , r 


t k = h^k, k = 0, 1, 2, S 


where h 1 and h 2 are positive integers. If h^, h 2 and the sampling time 
chosen so that h^T and h 2 T are small, then 


( r- 1 ) / 2 

Y(m) = E cos((r-2i)h. w-T) 

i =0 1 J 


( r -1 )/2 

= e cos ( ( 2i +1 )cj. h.T) 
i =0 J 1 


( r - 1 )/2 

= cos(u). h.T) [ e cos(2 0 ). i h,T)] 
J 1 i =0 J J- 


( r- 1 ) / 2 

- sin(cD.j h, T) [ e sin(2u . i h . T) ] 
J i=0 J 


w k T) 

( A5) 

(A6) 

(A7) 

(A8) 

T are 
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(r-l)h,T/2 

« cos(a). h.T) / cos(2w -t)/(h^T) dt 

j i o J 

(r-l)h.T/2 

- sin(o). h.T) / si n(2w -t)/(h^T) dt 

J 1 n ** 


= cos (o)j h^T) si n( ( r-1 ) W j h 1 T)/(2 h^) 

+ sin(<jjj h^T) [cos((r-l)wj h^T) - 11/(2 wj h^T) 

= (sin(r Wj h^T) - sin(cOj h^T ))/(2 uj h^T) 

| Y (m) | < 1 / (o)j hjT) 


< oo 


Similary, from Eq. (A5), it is shown that 


r Ji 

F 1 = Z (sin(j i cdjT) sin(j\ <d k T) + sin(j i w k T) si n( ( j p -j i )wjT) ) e 


= sin ( r to - h.T) z sin(u. i h.T) cosfo. i h.T) e 
J i i =0 K i J 1 


W 1 h i T 


r 

+ (1 - cos(r u. h.T)) Z sin (w. i h.T) sin ( u . i h.T) e 
J 1 i =0 k 1 j i 


(o j" a k )i 


<j.-o. r h.T r 

< e J K 1 {sin(r u. h..T) z [sin( (w. -a). )i h^t) + si n( • ) 

J i .j = Q ^ J 1 J 


r 

+ (l-cos(r uj hjT)) z [cos((a) k -ojj )i hjT) - cos( (w k +o>j )i hjT) ] } 


« e 


|o --a k | r h.T 
J fsinfr u>. 


{sin(ra)j h^T) [ (l-cos( (w k -a)j )r h^T) )/(2(a) k -wj )h^) 


(A9) 


(A10) 


(a r a k )T 


i h x T) ] 
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+ ( 1-cos ( (oj^+ojj) r h jT)/ (2(a)| < + a)j) h^T) ] 

+ (l-cos(roOj h^T)) [sin( (aj^-WjOr h 1 T)/( h^T) 


- s i n ( (ti> k +a3 j ) r h 1 T)/(2(<jd k +u) j .) h^J)]} 


l a r a kl r h l T 

|F 1 1 < ( 2/ (w^-Wj) + 2/ (a) k +u)j)) e /h^ 

For a finite data length, r h-j^T = constant, if T+0, then r+®. From 
(A10) and (A12), one obtains 


1 im 
T>0 
k*j 


This is also true if r is even. 


1 im 
T*0 
k*j 

Substitution of Eqs. (A13) and (A14).into (A5) yields 


- 0 

Similarly, it can be proved that 
F 2 /(2 x 2j .) = 0 


1 im 
TVO 
k* j 


V 

aj v lk 


= 0 


For Eq. (A4), similar procedures can be used to verify 


1 im 
T+0 

V ai V 2k = H ” 
aj T*0 

vf. V,, = lim 
bj lk ^ 

k*j 

k*j 

k*j 


1 im 

V i V k 

= 1 im W. wf 

T>0 

J k 

T+ 0 K J 

k*j 


k*j 




= 0 


(All) 

(A12) 

Eqs. (A6) , 
(A13) 

(A14) 

(A15) 

(A16) 


and 


(A17) 
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where 0 2 is a 2x2 zero matrix. Substution of Eq. (A17) into (A3) derives Eq. 
(Al). Note that a better approximation of Eq. (8) can be achieved if smaller 
T are used for a finite data length. 
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APPENDIX II 

For the controllability matrix W and the observability matrix V shown in 
Eqs. (15) and (16), the corresponding orthogonal matrices, W # and V^, can be 
derived as follows. First, observe that 

V # v = W W # = I 2 . (Bl) 

From Eqs. (17) and (25), it is shown that 

p 

Vg(r) Vj ( r) = E^ [sin 2 (j. w T) ( 1/X 1 ( r) + l/\ 2 (r)) + si n( j. u T)si n( ( j r -j . ) W T) 

(l/\ 1 (r) - l/\ 2 (r))]/2 = (l/2\ 1 (r)) e [sin 2 (j. u T) 

i=0 

r o 

+ sin (j i u)T)sin((j r -j.) w T)]+ 1/ ( 2x 2 ( r) ) e [sin (j^T) 


- sin (j.wT) sin((j r -j.)(oT)] (B2) 

If r is even, with the aid of Eqs. (21) and (29)- (31), one obtains 

r 2 

E^ [sin ( j i ajT) + sin(j iU T) sin( (j f .-j 1 )wT)] 

[sin(j iW T) + s i n ( ( j r - j.)uT)] 2 + 2 sin 2 (j p/2 U T) 


i=u 


r/2-1 

E 

1-0 


= E [4 cos 2 ( ( j . j/2) U T) sin 2 (j>T/2) + 2 sin 2 (jJ/2) 
i =0 1 r r r 

o r/2-1 , 

= 2 si n^ (j uT/2) [1 + E 2 cos*((j. - j /2) W T)] 

i=0 
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• 2 


r/2-1 


= 2 sin (j toT/2) [1 + E (1 + cos ((j - 2 j.) w T))] 

r i=0 r 1 


« 2 sin 2 (j r o)T/2) Xj (r) 


(B3) 


£ [sin 2 (j.coT) - sin(j.uT) s i n ( ( j - j.)coT)l 
i=0 1 r 1 


r/2-1 

S [sin(j.a)T) - si n( ( j - j.) w T)]‘ 
i =0 1 r 1 


r/2-1 ? 

E [4 cos^(j u T/2) sin^( (j . - j /2) W T)] 
i=0 r 1 r 


.2 


r/2-1 


2 cos (j coT/2) E [1 - cos((j - 2j.) w T)l 
i=0 r i 


- 2 cos 2 (j r o>T/2) \ 2 (r) 


(B4) 


If r is odd, with the aid of Eqs. (21) and (29)-(31), one arrives 
^ 9 

E [sin- (j.coT) + sin(j,(oT) sin((j - j.) w T)l 
i =0 1 1 r l 

( r- 1 ) / 2 

= £ [sin(j.toT) + sin((j - j.)a)T)]^ 

i=0 1 r i 


2 ( r- 1 ) / 2 

= 2 sin (j coT/2) e [1 + cos((j - 2j.)u,T)] 

i =0 


= 2 sin 2 (j r u>T/2) ^(r) 


(B5) 


. 2 


E [sin (j.coT) - sin(j. w T) sin ((j r - j.) ojT) ] 
i=0 1 1 r i J 
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( r- 1 ) / 2 2 

= E [sin(j.a)T) - sin ( ( j - j , ) ojT ) ] 
i=0 1 r 1 

o (r-l)/2 

= 2 cos* ( j toT/2 ) z [l-cos((j -2j.)a)T)] 

i=0 

= 2 cos 2 (j r o)T/2) \ 2 (r) (B6) 

Substitution of Eqs. (B3)-(B6) into (B2) yields 

vj(r) V^r) = [2 sin 2 (j^T/2) ^ (r)]/^ (r) + [2 cos 2 (j r a,T/2)\ 2 (r)]/2\ 2 (r) 

= 1 (B7 ) 

Simil arly 

Vj(r) V 2 (r) = 1 (B8) 

Next, from Eqs. (18) and (25), it is shown that 


Vg(r) V 2 (r) = - e^ [sin(j.a)T) cos(j.wT) + sin((j p - j i )a>T)cos(j.< 1 )T)]/2\ 1 (r) 


+ E [- sin(j.u>T) cos(j i u)T) + sin((j - i . )<*>T) cos(j .u>T)]/2\ ? (r) 
i=n L 

U (B9) 


If r is even, with the aid of Eqs. (21) and ( 29 ) - ( 3 1 ) , one obtains 


E [sin(j.a)T) cos(j.wT) + s i n ( ( j - j, )u>T) cos(j. w T)] 
i =0 


r/2-1 

E [sin(j.uT) + sin( (j - j . )ojT) ] [cos( j - ojT) + cos((j r - j,)o>T)l 

i =0 r i 
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+ 2 sin ( j r a)T/2 ) cos(j r wT/2) 


r/ 2-1 

= E 4 si n( j mT/ 2) cos((j.-j /2) U T) cos(j w T/2) cos((j.-j /2) U T) 
i=0 r 1 r r r 


+ sin(j r uT) 


r/2-1 7 

= s i n ( j gjT ) [1 + j 2 cos ((j. - j /2) W T)] 
r i =0 i r 


r/2-1 

= sin(j r uT) [1 + E (1 + cos(’(j r - 2 j^wT)] 


= sin(j rW T) x^r) 


r 

E [- sin(j.wT) cos(j. w T) + sin ((j - j. ) w T)cos( j .«T)] 

i =0 1 1 r 1 1 


r/2-1 

= Z [sin( (j - j . )o)T) - sin(j. w T)] [cos(j. u T) - cos((j - j,)u>T)] 
i=0 


r/2-1 

= E 4 [sin( (j/2 - j. )uT)cos(j mT/2)] [- sin( j r u,T/2)sin( (j . - j/2) 
i =0 


r/2-1 ? 

= sin(j o)T) E 2 sin ((j /2 - j.) u T) 
r i=0 r 1 


r/2-1 

= sin(j U T) E [1 - cos ( ( j - 2 jDwT)] 
r i =0 r i 


= sin(j r<J jT) x 2 (r) 


(B10) 


toT)l 


(Bll) 
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Substitution of Eqs. (BIO) and (Bll) into (B13) yields 
vj(r) V 2 (r) = - [s1n(j r a>T)\ 1 (r)]/2\ 1 (r) + [sin(j r wT)\ 2 (r)]/2\ 2 (r) = 0 

( B1 2 ) 

This is also true if r is odd. Similarly, it can be proved that 

vJ(r)V 1 (r) = 0 (B13) 

Observation of Eqs. (15), (23), (B7), (B8), ( B 1 2 ) and (B13) leads to 

V # V = I 2 ( B14) 

Similar procedures can be used to verify 

WW # = I 2 (B15) 


Note that, if a = 0, from Eqs. (15), (16), (23), (24), it can be proved 

that 


(VV # ) T = VV # , (wV = W # W 


(B16) 


From Eqs. (B14) - (B16), it shows that and W # are the pseudo inverse 
matrices of V and W respecti vely. However, for a * 0, and are not the 
pseudo inverse matrices. This raises doubts about the uniqueness of the modal 
parameters identified. The effect of this non-unique parameters needs further 
study in the future. 
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Table 3.1 The Global Minimum of the 
Three Hump Camel Function 


No. 

X 1 

x 2 

fU}, X 2 ) 

1 

0.000 

0.000 

0.00 


Table 3.2 The Global Maxima of the Himmelblau Function 


No. 

X 1 

x 2 

F(x 1 ,x 2 ) 

1. 

3.0000 

2.0000 

0.00 

2. 

3.5844 

-1.8481 

0.00 

3. 

-3.7793 

-3.2832 

0.00 

4. 

-2.0851 

3.1313 

0.00 


Table 4.1. Percentage Error of Estimated Modal Frequency 


Damping 

X^xratios 

\ Noise 

Frequency \ 
(rad/sec) \ 

0.3% 

2% 

0% 

5% 

10% 

30% 

0% 

5% 

10% 

30% 

3 

0.86 

1.10 

1.33 

2.28 

0.75 

1.00 

1.26 

2.23 

7 

0.03 

0.33 

0.66 

1.95 

0.04 

0.42 

0.65 

2.03 

11 

0.18 

0.26 

0.59 

1.63 

0.17 

0.39 

0.82 

1.92 

17 

0.01 

0.17 

0.33 

0.99 

0.01 

0.25 

0.53 

1.72 

23 

0.02 

0.08 

0.13 

0.31 

0.01 

0.12 

0.21 

0.48 

31 

0.01 

0.17 

0.32 

0.84 

0.02 

0.57 

0.96 

1.60 

37 

0.01 

0.02 

0.06 

0.35 

0.02 

0.42 

0.93 

1.70 

43 

0.01 

0.18 

0.36 

0.93 

0.00* 

0.79 

1.17 

1.54 

53 

0.00* 

0.03 

0.06 

0.14 

0.01 

0.13 

0.21 

0.36 

67 

0.00* 

0.13 

0.27 

0.63 

0.00* 

0.90 

0.95 

1.01 


*after numerical truncation 


Table 4.2. Percentage Error of Estimated Modal Damping 


Damping 
\ ^\ratios 

\Noise 


0 . 

32 


22 

Frequency \ 
(rad/sec) \ 

0% 

52 

102 

302 

02 

52 

102 

302 

3 

128.20 

* 

★ 

* 

17.21 

28.29 

73.22 

★ 

7 

3.51 

it 

it 

* 

0.33 

53.06 

★ 

* 

11 

20.79 

★ 

★ 

* 

3.34 

25.72 

52.03 

* 

17 

0.21 

4.88 

11.03 

60.52 

0.13 

0.18 

0.46 

15.92 

23 

13.19 

1.93 

14.73 

65.56 

0.81 

7.42 

14.42 

35.66 

31 

10.53 

1.42 

17.92 

* 

2.68 

1.66 

13.07 

68.62 

37 

1.99 

24.73 

46.44 

126.50 

0.11 

13.03 

9.71 

27.50 

43 

0.13 

0.77 

7.07 

* 

0.89 

12.90 

34.60 

85.88 

53 

0.49 

15.97 

32.01 

91.41 

0.19 

14.54 

25.78 

55.24 

67 

1.65 

9.09 

5.23 

76.90 

0.31 

40.22 

58.45 

88.26 


* The estimated modal damping is 0 % 


cn 

cn 
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Table 4.3 Percentage Error of Estimated Modal Frequency 
of the MAST Beam Structure 


Damping 

ratio 

0% 

0.3% 

2% 

n. Noise 

Level 

0% 

5% 

10% 

0% 

5% 

10% 

0% 

5% 

10% 

Freq 

(Hz) \ 










1.4222 

0.00# 

0.01 

0.02 

0.01 

0.00# 

0.01 

0.00# 

0.00# 

0.01 

1.4222 

0.02 

0.06 

0.16 

0.03 

0.06 

0.15 

0.03 

0.09 

0.18 

8.5545 

0.21 

0.21 

0.20 

0.23 

0.23 

0.23 

0.38 

0.40 

0.42 

9.4954 

0.82 

0.80 

0.78 

0.84 

0.81 

0.79 

0.97 

0.95 

0.90 

9.4954 

0.29 

0.30 

0.30 

0.30 

0.31 

0.31 

0.43 

0.43 

0.42 


# after numerical truncation 


Table 4.4 Percentage Error of Estimated Modal Damping 
of the MAST Beam Structure 


Damping 

ratio 


0% 



0.3% 



2% 


Noise 
n. Level 

Freq n. 

(Hz) \ 

0% 

5% 

10% 

0% 

5% 

10% 

0% 

5% 

10% 

1.4222 

★ 

★ 

★ 

0.33 

3.33 

6.67 

0.00# 

0.50 

1.00 

1.4222 

★ 

* 

★ 

6.0 

10.00 

13.33 

0.00# 

0.50 

1.00 

8.5545 

★ 

★ 

★ 

★ 

* 

★ 

41.0 

41.50 

42.00 

9.4954 

★ 

★ 

★ 

★ 

* 

★ 

75.0 

77.00 

78.50 

9.4954 

★ 

★ 

★ 

83.3 

96.67 

★ 

38.5 

42.00 

44.00 


* The estimated modal damping is 0% 

# after numerical truncation 



Fig. 3.1. The Steps of the Optimization Algorithm 
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Fig. 3.2 


A Surface Plot of the Three Mump Camel function 
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3.4 A Surface Plot of 


the llimnel blau function 



Fig. 3.5 Contours of the Surface of the Himmelblau Function 



Fig. 4.1 Simulated Free Impulse Response Data 


rv> 


Amp! i tude 


Noise Level : 30? 


UD . 





Time 


f i o . 4.2 Simulated Measurement White Noise 
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Fig. 4.3 Cost as a Function of Updated Frequencies 
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Fig. 4.4. The MAST truss beam structure 
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Fig. 4.5 Mode shape of the first mode of the MAST truss beam structure. The 
solid bars indicate the displacement in the x-direction. The empty 
bars indicate the displacements in the y-direction. There are two 
measurements on each bay for each direction. 
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Fig. 4.6 Mode shape of the second mode of the MAST truss beam structure. 

The solid bars indicate the displacements in the x-direction. The 
empty bars indicate the displacements in the y-direction. There 
are two measurements on each bay for each direction. 
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Fig. 4.7. The simulated free impulse response data with 10% noise and zero damping 
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Fig. 4.8 


The cost function J for 
their frequencies with a 


the first mode as the 
fixed zero damping for 


projection filters update 
noise free case. 
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